Eccentricity Evolution of Extrasolar Planets 



Eccentricity Evolution of Extrasolar Multiple Planetary Systems 
due to the Depletion of Nascent Protostellar Disks 

M. Nagasawa^'^ D. N. C. Lin\ and S. Ida^ 

mnagasawaSmail . arc . nasa . gov , linSucolick . org , idaOgeo . titech .ac.jp 

ABSTRACT 

Most extrasolar planets are observed to have eccentricities much larger than 
those in the solar system. Some of these planets have sibling planets, with com- 
parable masses, orbiting around the same host stars. In these multiple planetary 
systems, eccentricity is modulated by the planets' mutual secular interaction as 
a consequence of angular momentum exchange between them. For mature plan- 
ets, the eigenfrequencies of this modulation are determined by their mass and 
semi-major axis ratios. But, prior to the disk depletion, self gravity of the plan- 
ets' nascent disks dominates the precession eigenfrequencies. We examine here 
the initial evolution of young planets' eccentricity due to the apsidal libration 
or circulation induced by both the secular interaction between them and the self 
gravity of their nascent disks. We show that as the latter effect declines adiabati- 
cally with disk depletion, the modulation amplitude of the planets' relative phase 
of periapse is approximately invariant despite the time-asymmetrical exchange of 
angular momentum between planets. However, as the young planets' orbits pass 
through a state of secular resonance, their mean eccentricities undergo systematic 
quantitative changes. For applications, we analyze the eccentricity evolution of 
planets around Upsilon Andromedae and HD 168443 during the epoch of proto- 
stellar disk depletion. We find that the disk depletion can change the planets' 
eccentricity ratio. However, the relatively large amplitude of the planets' eccen- 
tricity cannot be excited if all the planets had small initial eccentricities. 
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1. INTRODUCTION 

Many extrasolar planetary systems have been discovered recently using the radial ve- 
locity technique (Marcy & Butler 2000). The basic assumption is that the spectroscopic 
variations observed in some targeted stars are due to the Doppler shift associated with the 
reflex motion of stars with unseen companions. The mass of these companions depends on 
the poorly known inclination of these systems. But, unless the orbits of these systems are 
highly inclined (Stcpinski & Black 2001), the inferred masses of the companions arc compa- 
rable to that of Jupiter. The observed spectra of some targeted stars indicate the presence of 
multiple companions around them. For example, three companions are found to orbit around 
Upsilon Andromedae. The long term stability of this system requires their inclination to be 
sufficiently small such that the masses of these companions are no more than a few Jupiter 
mass (Mj) and much smaller than that of their host star (Laughlin & Adams 1999; Rivera 
& Lissauer 2000; Stepinski, Malhotra, & Black 2000; Ito & Miyama 2001; Lissauer & Rivera 
2001). Two companions are also found around HD168443 with minimum masses which are 
an order of magnitude larger than Mj. Although the masses of these companions would be 
a modest fraction of a solar mass (M©) if this system is viewed nearly face on, such com- 
pact hierarchical stellar systems have not be seen before. Thus, we follow the conventional 
practice to refer these multiple companions as planets. 

In the limit that the planets' masses are substantially smaller than that of their host 
star, their mutual secular perturbation induces them to exchange angular momentum while 
preserving their energy. (Planets in mean motion resonances also exchange energy on com- 
parable time scales). Around Ups And and HD168443, the planets' orbits are not in mean 
motion resonances so that their semi-major axes are conserved while their eccentricity and 
longitude of periapse modulate over some characteristic secular time scale (Murray & Der- 
mott 1999). 

But, the secular perturbation between the planets has not always been sustained at the 
present level. During the epoch of their formation, protoplanets are embedded in protostellar 
disks which have been found around most young stellar objects (cf. Haisch, Lada, & Lada 
2001). The mass and temperature distribution of these disks are very similar to those inferred 
from the minimum mass nebula model for the solar system (Beckwith 1999). The self gravity 
of these disks can induce the orbits of planets formed within them to precess at a rate faster 
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than that due their mutual perturbation. Consequently, the rate of angular momentum 
transfer between the interacting planets is suppressed. 

In the solar system, the initial contribution of the solar nebula to the total gravita- 
tional potential dominates the precession frequency of asteroids and comets over that due 
to the secular perturbation induced on them by the giant planets. The planets also undergo 
precession induced by the disk gravity and other planets' secular perturbation. In general 
the precession frequencies of these celestial bodies do not equal each other. But as gas is 
depleted in the solar nebula, its self gravity weakens. The total precession frequencies of 
both the asteroids and the planets declines, though at a different rate. When the precession 
frequency of asteroids with some semi major axes coincides with that of a major planet, they 
enter into a state of secular resonance. In this resonance, the eccentricity of the asteroids 
is either excited or damped monotonically, depending on their relative longitude of perias- 
tron passage with respect to that of the planet. As this secular resonance sweeps across the 
solar system, large eccentricity may be excited among some small celestial bodies (Ward, 
Colombo, & Frankhn 1976; Heppenheimer 1980; Nagasawa, Tanaka, & Ida 2000; Nagasawa 
& Ida 2000). 

In this paper, we examine the effects of disk depletion on the eccentricity evolution of 
the planetary systems around Ups And and HD 168443. Through such an investigation, we 
hope to infer the kinematic properties which these planets are born with and thereby cast 
constraints on their formation process. In §2, we briefly discuss the precession due to the 
secular interaction between planets and that due to the self gravity of the protostellar disk. 
The conditions for secular resonance are discussed. We introduce, in §3, a working model 
and describe the method we used to analyze the eccentricity evolution. In §4, we present the 
results of some calculations. We use these numerical results to demonstrate the evolution of 
these systems with Hamiltonian contour maps. Based on these results, we infer, in §5, some 
implications on the planets' orbital eccentricity shortly after their formation and while they 
are still embedded in their nascent disks. Finally, we summarize our results in §6. 

2. SECULAR INTERACTION IN MULTIPLE EXTRASOLAR 

PLANETARY SYSTEMS 

2.1. Current Orbital Properties 

In this paper, we focus our discussion on the planets around Upsilon Andromeda (c 
and d) and those around HD168443. The decomposition of Upsilon Andromeda's spectra 
indicate that there are three planets orbiting around it (Butler et al. 1999). Inferred orbital 
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elements^ of the Upsilon Andromeda planets are shown in Table 1(a). The semi-major axis, 
the eccentricity, the longitude of pcriastron, mass of the planet, and periastron passage time 
(JD) are denoted by a, e, zu, M, and Tperi, respectively. The subscripts b, c, and d represent 
the values for the individual planets. Throughout this analysis, we assume that these systems 
are viewed edge-on and their masses correspond to their minimum values. Provided their 
orbits are coplanar, our analysis is independent of the planets' inclination. Our approach 
becomes inadequate for the limiting cases of nearly face-on orbits where the masses of the 
companions are comparable to their host stars. 

Around Ups And, the eccentricity of the innermost planet b is essentially undetectable. 
The orbit of planet b is likely to be circularized during the main sequence life span of the host 
star Ups And by the tidal dissipation within the planet's interior as expected for Jupiter-like 
extrasolar planets with semi- major axis less than 0.05 AU (Rasio et al. 1996). With its low 
mass and small semi-major axis, planet b does not contribute significantly to the dynamical 
evolution of the system (Mardling & Lin 2003). Although the outer planets (c and d) still 
exert secular perturbation on planet b, the cumulative effect on its eccentricity modulation 
is limited. Under the present configuration, the secular interaction of planet b with planets c 
and d is weakened by the rapid precession due to the post-Newtonian relativistic correction 
in the gravitational potential of the host star (Mardling & Lin 2003). Neglecting the secular 
perturbation due to planet b, the orbital evolution of planets c and d obtained from the direct 
orbital integration of the full equation of motion is shown in Figure la. In the numerical 
integration, the calculation is started with the eccentricity ratio x = ec/e^ and the relative 
longitude of periastrons between planets c and d, r) = zUc — vjd, based on their observed 
values 0.66 and ~ —0.06 radian, respectively. We compute the evolution of the system for 
over 2 x 10^ years. The equi-Hamiltonian contours (see Appendix D) are also shown. For 
planets c and d around Ups And, the massive planet d has larger eccentricity than smaller 
planet c. The longitudes of pcriastron of planet c and d are always close to each other. When 
the planets arc in librating (closed) track in (x-f]) diagram, the planetary system tends to be 
stable for a long time. The stabilities of this system around Ups And are well investigated 
(Laughlin & Adams 1999; Rivera & Lissauer 2000; Stepinski et al. 2000; Barnes & Quinn 
2001; Ito & Miyama 2001; Lissauer & Rivera 2001; Chiang, Tabachnik, & Tremaine 2002; 
Mardhng & Lin 2003). 

Around HD168443, two massive planets are inferred from the radial velocity curves 
(Marcy et al. 2001). The best-fit orbital elements are shown in Table 1(b). We also numeri- 
cally integrate the present-day orbital evolution of the planets b and c around HD168443. In 
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contrast to the planetary system around Ups And, the lower mass planet b around HD 168443 

system, has an eccentricity which is more than twice that of the massive planet c. The orbits 
are integrated over 2.5 x 10^ years, starting with x = 2.65 and rj ~ 1.92 radian. These planets 
evolve in circulating (open) track in the [x-f]) diagram (Fig. lb) and t] evolves from — tt to vr 
over a period of about 1.6 x lO'' years. This system too is stable despite the large magnitude 
of eccentricity of planet b. The origin of such large mass and eccentricity of planets in this 
system has not been addressed previously. 



2.2. Secular Perturbation between Two Planets 

The secular interaction induces eccentricity modulation between planets c and d around 
Ups And (Laughlin & Adams 1999; Rivera & Lissauer 2000; Ito & Miyama 2001; Mardling 
& Lin 2003). For presentation purpose, it is useful to briefly recapitulate the analysis of 
secular interaction between two planets. It is customary to consider the secular (long-term) 
evolution of the interacting planets' orbits using a disturbing function (Murray & Dermott 
1999). To the lowest order, the modulation of the eccentricity {ec,d) and longitude of periapse 
(zucd) of some planets c and d can be approximated by 

= lc,dCed,cSm{wc,d - ^d,c) + K,d: (1) 

^ = led {l + C cos(.z7,,, - + + ^c,d, (2) 

where r = t/tc, t is the real time, t,, = {4/nf.){M^,/Md){a(i/ai,y /b^y^, bi^\a) are the Laplace 
coefficients with a = ac/aa, M^, M^, and are the mass of planets c, d, and the host 
star respectively, Oc and aa are the semi-major axes of planets c and d respectively, nc,d = 
{GM^/al^^y^'^ is the mean motion of planet c or d, ^c,d = {o-d o.c,d)^^'^{Mc/ M^^d) so that 
7c = 1, C = — ^3/2/^3/2 (Mardhng & Lin 2002), A^^^ is the change rate of eccentricity 
due to the non axisymmetric torque, A^^ is the change rate of vo due to the non axisym- 
metric torque, and Nc^ is axisymmetric modification of the gravitational potential. When 
flc/flrf -C 1, tc = (4/3nc)(M^,/M(;)(a(i/ac)^ and C = —bad'^ad- Frequently used variables are 
tabulated in Table 2. The planets' secular interaction induces both axisymmetric and non 
axisymmetric perturbations on each other. While only the non axisymmetric component of 
the planets' secular perturbation on each other can lead to angular momentum exchange 
and ec,d modulations in equation (1), both axisymmetric and non axisymmetric perturba- 
tions contribute to the evolution of va^d through the first and second term on the right hand 
side of equation (2) respectively. These linearized equations of secular motion, though only 
accurate to the first order, are useful for obtaining analytic-approximation solutions which 
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highlight the dominant physical effects. For the numerical calculation of the planets' orbits 
(see §4), we use the exact equation of planets' motion. 

Planet-disk interaction also leads to both axisymmetric modification of the gravitational 
potential (A^cd) and non axisymmetric torque (A^'J) (Goldreich & Tremaine 1980, 1982; Lin 
& Papaloizou 1986 a, b, 1993). For the latter effect, planets excite waves in the disk which 
carry angular momentum. The dissipation of these waves, anywhere in the disk, would 
result in a finite torque (Papaloizou & Lin 1984). In principle, the net torque vanishes in 
the invisid limit. However, their amplitude grows and steepens into nonlinear shocks as the 
waves propagate away from the location where they are launched, leading to an effective 
torque (Savonije, Papaloizou, & Lin 1994). The evolution of the axisymmetric potential 
modifies the precession of w but do not directly infiuence the eccentricity of the orbits (see 
eq. [2]). While the non axisymmetric torque may induce monotonic changes in e (Chiang 
& Murray 2002) and w (see eq. [1]), the magnitude and sign of A^^ depends sensitively 
on the disk structure. Interaction between the embedded planets with disk gas through 
corotation resonances damps the eccentricity while that through Lindblad resonances excites 
the eccentricity (Goldreich & Tremaine 1980). For protoplanets with mass less than a few 
times that of Jupiter and modest eccentricity, gas may flow in the vicinity of their orbits 
such that the eccentricity damping effect of the corotation resonances is stronger than the 
excitation effect of the Lindblad resonances (Goldreich h Tremaine 1980). But protoplanets 
with masses an order of magnitude larger than that of Jupiter may open relatively wide gaps 
in protostellar disks. In this limit, the protoplanets' corotation resonances may be cleared of 
disk gas such that their eccentricity may be excited (Artymowicz 1993; Papaloizou, Nelson, 
& Masset 2001; Goldreich & Sari 2002). For precession, the contribution from the non 
axisymmetric torque may be weaker than that due to the disk's axisymmetric contribution 
to the total gravitational potential in the limit that the torque resulting from the planets' 
interaction with the interior and exterior regions of the disk are balanced or in low- viscosity 
and thin disks where planets with relatively low masses can open wide gaps. However, the 
non axisymmetric torque may lead to significant contributions over time in relatively massive 
disks. 

In order to analyze the contribution of each effect, we adopt a piece meal approach. In 
this paper, we focus our attention on the evolution of planetary orbits due to the changes in 
the axisymmetric disk potential. For the departure from a point-mass potential, the apsidal 
motion of planets c and d are included in Nc,d (see Appendix A). In this first step, we 
neglect the effects due to the non axisymmetric torque by setting A^'J = in equations (1) 
and (2). The extension of our discussion to the limit of finite non axisymmetric torque will 
be presented in a future contribution. 
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2.3. Precession Frequencies 

To the lowest order, the solution of equations (1) and (2) can be expressed as 

ec,d exp(iroc,d) = Ac4exp{i{git + + Bc,dexp{i{g2t + P2)}, (3) 

where Ac^d and fic.d ai'e the oscillation amplitudes, (3i and (32 arc the phase angles (e.g., 
Brouwer & Clemence 1961). The individual planet's longitudes of periastrons precess with 
two independent eigenfrequencies, 

51,2 ^^[gc + gd± [{gc - gaf + ^gla] ^'^ } , (4) 

where 

gc,d = — -, gcd = —f— [McMd) ' . (5) 

Prom equation (2), we find that gc would be the precession frequency of planet c if it is 
massless and is subjected to the secular perturbation of a massive planet d with a circular 
orbit. Similarly, would be the precession frequency of planet d if it is massless and is 
subjected to the secular perturbation of a massive planet c with a circular orbit. 

In the limit that A^^d a-nd 3^,^, have comparable magnitudes and the two eigenfrequen- 
cies, g\^2 have sufficiently large differences, the relative longitude of the two planets would 
circulate. Precessional degeneracy occurs when g\ = g2 so that the two planets precess 
synchronously regardless of the amplitude of their e and w. In this degenerate state, the 
relative phase of the two planets' longitude of periapse passage, r] = Wc — ^d, retains a 
constant value and the two planets are in a state of secular resonance. In this secular reso- 
nance, angular momentum is monotonically transferred from one planet to another because 
the relative orientation of their periastrons is maintained indefinitely. To lowest order, the 
planets' individual energy is conserved. With the exception of the special configuration in 
which = or ±7r, this monotonic transfer generally leads to eccentricity excitation of 
the angular momentum losing bodies and eccentricity damping of the angular momentum 
gaining bodies. In celestial mechanics, secular resonance is thought to be important in the 
current eccentricity distribution of the asteroids (Williams 1969). 

Prom equation (4) , we find that the magnitude of the difference between gi and g2 is 

1 /9 

^5-12 = \gi - g2\ = [{gc - gaf + ^gcd\ = 

which vanishes on the exact center of secular resonance (where gi — g2) only if gc — gd and 
gcd — 0. The former requirement corresponds to a necessary resonance condition = 



1 - 7d + AiV 



4^, 



cd 



(6) 
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where 

A = l-7rf + AAr (7) 

is the precession rate induced by axisymmetric component of the perturbed potential (see 
eq. [9] below). The latter requirement is satisfied if either Mc = or — (see eq. 
[5]). Equation (6) implies that planets with comparable masses cannot have precessional 
degeneracy with gi — g2 (Kinoshita & Nakai 2000). 

The present value of 7^ is 0.29 for planets c and d around Ups And. Today, in the 
absence of any residual disk (A^c,<i — 0), 92 9d < 9i 9c ^-nd the outer two planets of Ups 

And are not in a state of precessional degeneracy. Nevertheless, their present orbits can be 
approximated by equation (3) with the magnitude of A^^d much smaller than that of -Bed- 
Thus, the two planets primarily precess with eigenfrequency gi and their relative longitudes 
of periastrons librate over a restricted range of phases. This phase lock is equivalent of a 
state of secular resonance. This resonant interaction is the result of a dynamical feedback 
through the planets' secular interaction. 



2.4. Three Classes of Relative Orbits 

In order to illustrate the importance of this feedback effect, we substitute x = Cc/cd and 
T) = zUc — zud, equations (1) and (2) reduce to 

dx 

— = C(l + ^dX^)smri, (8) 
^ = (1 - 7rf) + iC/x){l - -fdX^)cosr] + AN = Ai + A2ix)cosr], (9) 

(IT 

where AN = Nc — Nd, Ai is given in equation (7), and A2{x) = {C/x){l — ■jdx'^) is introduce 
for notational simplicity. There are three families of relative orbits which can be illustrated 
below with the linearized approximation solutions of equations (8) and (9) (see Appendix 
B). We show below that the relative magnitude of Ai and A2 determine the nature of the 
orbits. 

1) CirculaiAon. In these solutions, the magnitude of x modulates about some values Xq 
such that X = xq + 6x{t) with an amplitude |5a;(r)| ^ xq for all values of t] which ranges 
between and 27r. In the limit that \Ai\ ^ |^2|j f] decreases monotonically (because Ai < 0) 
while 5x oscillates (see solutions in eqs. [B7] and [B8] in Appendix B.2). To the lowest order 
these solutions reduce to 

C 

77 ~ AiT, 5x ~ — COS^iT. (10) 
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Near the secular resonance where \Ai\ is relatively small, the amplitude of 5x becomes large. 

2) Libration. In the opposite limit that l^il -C IA2I, the non axisymmetric secular 
interaction between the planets is important. There are stationary points in the (x-r)) plane 
which center on the values oi x — Xm and 77 = or tt (see Appendix B.l). Around these 
points, there are orbits with both small amplitude modulation such that 

e = = eosma;r, 77 = 77oCOsa;T, (llj 

where 

C 

r)o = ±60 and cu = ± — (l + 7^0;^) (12) 

is the oscillation frequency (see Appendix B.l). The dimensionless amplitude of these libra- 
tional orbits eo <^ 1. For these orbits, although gi ^ g2, secular interaction induces them 
to process at similar frequencies that their relative longitude of periapse passage is always 
approximately aligned or anti aligned. Thus, these planets are effectively in a state of secular 
resonance. Note that because the libration is centered around 77 = or tt, there is no effective 
angular momentum transfer despite the phase lock. 

3) Excitation. For systems with \ Ai\ < IA2I, there are also orbits in which x and 1] have 
very different values as Xm and (or tt) respectively. In these cases, i] would evolve rapidly 
to a phase angle 

such that drj/dr is reduced to zero (see eq. [9]) and the two planets become phase locked. 

For all non zero (or tt) values of r]i , the monotonic increases/decreases of x correspond 
to eccentricity excitation/damping, analogous to the situation of precessional degeneracy. In 
this state of near secular resonances, angular momentum is monotonically transferred from 
one planet to the other, resulting in a monotonic evolution of x. The modification of x 
in turn leads to an evolution in r). Along the path of x and 771 evolution, gi ^ g2 and the 
planets would libratc about their evolving guiding center. For the special cases when the two 
planets enter the resonance, the second order solution reduces to that in equation (11) and 
all orbits become instantaneously librational, even though eo and 770 may become arbitrarily 
large. At the center of resonance where Ai, A2, and A1/A2 all vanish, 7/1 = n/2 and Xm 
evolves exponentially. We discuss the evolution of these systems in Appendix C. 
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3. MODELS 

The above analytic approximation is useful for isolating the three family of orbits along 
with two stationary points in the {x-r]) diagram. These orbits generally follows the contours 
of equi Hamiltonian map (see Appendix D). However, the topological evolution of the 
Hamiltonian map alone is insufficient for the determination of the planets' orbits during 
the depletion of the disk. We carry out, below, numerical integration of the full equation 
of motion for the planets subject to the potential of their host star and nascent disks. In 
this section, we briefly describe a model prescription with which we examine the passage of 
librational degeneracy during the epoch of disk depletion. 

3.1. Planetary Formation Scenarios and Disk Model 

We first discuss the physical process of disk depletion. According to conventional theo- 
ries, planets are formed through the condensation of grains which grow to planetesimals via 
cohesive collisions (Hayashi, Nakazawa, & Nakagawa 1985; Lissauer 1987; Wetherill 1990). 
Upon attaining a sufficiently large mass, planetesimals accrete gas (Mizuno 1980; Boden- 
heimer & Pollack 1986; Pollack et al. 1996). Eventually, their growth is terminated when 
protoplanets can tidally induce the formation of a gap near their orbit (Goldreich & Tremaine 
1980; Lin & Papaloizou 1980, 1993; Takeuchi, Miyama, & Lin 1996). 

Thereafter, gas in the inner of the disk continues to diffuse inward as it loses angular 
momentum to the planet. Since gas replenishment is cut off by the formation of the gap, the 
inner region is depleted well before the outer region of the disk. For the present discussions, 
we assume that the mass in the inner regions of the disk becomes negligible by the stage 
when the second planet is fully grown. Thus, in our analysis of planets' dynamical evolution, 
we only need to consider their interaction with the disk region extended outside the outer 
planet. Both the gravitational and tidal influences of the disk on the outer planet are much 
more intense than on the inner planet because the former is much closer to where most of 
the gas is distributed. 

When the planets' orbits lie in the same plane as their nascent disks and their distance 
to the nearest edge of their nascent disks is greater than the disks' scale height, a thin- 
disk approximation is sufficient for the computation of the disks' gravitational potential. 
Following the approach of Ward (1981), the self-gravitating potential for a disk model with 
a surface density profile E = Eo(ro/r)'^ (with Eq and Tq being some fiducial values) can be 
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expressed as 



V{r) = 27rG'Eor 




In the above expression, r is the instantaneous location (rather than the semi major axis) 
of the planet. We also assume that the disk has an inner edge which is located at Tedge, 
interior to which E = 0. The coefficient An = {(2n)!/2^"(n!)^}^. For illustration purpose, 
we consider the minimum mass solar nebula model (Hayashi 1981), in which k = 3/2. The 
general applicability of this phenomological model for protostellar disks around different 
stars is questionable. In the absence of a reliable prescription, we parameterize disk models 
with a range of values in k (see §4.3). We assume the disk is extended to infinite and the 
semi-major axes of planets do not change as a consequence of disk depletion or planets' 
secular interaction with each other. The results to be presented below show that while the 
overall dynamical evolution of the planets' orbit is insensitive to the detailed disk model, the 
planets' extrapolated initial eccentricity ratio does depend on the mass of the disk gas near 
their orbits. 



3.2. Planetary System and Prescription for Disk Depletion 

For computational simplicity, we consider a planetary system consists of two planets 
with coplanar orbits. In the case of the Upsilon Andromeda system, we are concerned with 

gravitational interaction between planets c and d and neglect the contribution by planet 
b. We also neglect the influence of the planet's tidal interaction with the disk. Only the 
axisymmetric component of the disk potential is taken into account. Mutual inclinations 
between the two planets and between the planet and the disk arc neglected. The inclination 
of the orbital plane is denoted by i. We take surface density of the disk as 5 x sin i times that 
of the minimum mass solar nebula model. Only the ratios, rather than the individual values, 
of the planets' and the disks' masses determine the planets' motion. Thus, our results do 
not depend on the inclination of the system (except for small changes in the mean motion 
and the long term dynamical stability for exceptional small values of i). 

For models with different values of k, we scale E such that its value at the disks' inner 
edge Tedge remains invariant. As we have indicated above, inner region of the disk is likely to 
be severely depleted prior to the emergence of the second planet. Thus, we assume the disk 
depletion to proceed in an inside out manner by adopting a model in which the location of 
nebula inner edge ("Tcdgc) is prescribed to expand outwards as 'rcdgc(^) — ''"edge 

(t = 0)+rit/tAN 

with Ti = lAU while Eq is kept constant. In order to demonstrate that the passage of the 
secular resonance does not depend on the detailed prescription of the depletion process or 
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the structure of the disk, we also consider a model in which the surface density of the disk 

is assumed to decline by an identical reduction factor everywhere such that Eo(t) = Eo(t = 
0) cxp(— t/tAJv) with Tedge =constant. In the limit that a^/a^ ^ 1 and ac^/r^dge ^ 1, we 
find that the two prescriptions of the disk depletion do not lead to a significant difference in 
the numerical results because AA^ oc S/r^^g^. For the approximate prescription of AA^ (see 
eq. [D2]), the magnitude of Ai = 1 + AA^ — changes on a time scale ta^ ~ ^aa? and the 
condition of the secular resonance = is satisfied when 

^orf ^,5{j^-l)M, f^_n^\-\ ^^^^ 



edge " ^ 

This condition can be attained with both prescriptions. 

The orbital evolution of the planets can be directly integrated numerically, including 
the point-mass potential of the host stars and sibling planets as well as the contribution from 
the disk potential where the gravity can be expressed as 

dr 

where A; = 3/2 is assumed. 
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4. RESULT OF NUMERICAL COMPUTATION 

4.1. Evolution of Librating Planets Around Ups And 

Using the equi-Hamiltonian contour maps, we illustrate the orbital evolution of planets c 
and d during disk depletion. Figure 2 shows the time evolution of the equi-Hamiltonian map 
of Upsilon Andromeda system. In these calculations, we adopt the mass ratio Mc/Md = 0.502 
based on the assumption that they occupy the same orbital plane. The nebula edge retreats 
from inside to outside (from panel 1 to panel 5). When the nebula edge is at 4.2AU (panel 
1), the orbits librate around {x = 10, rj = 0) or {x = 0.4, r] = ±7r). The precession speed of 
the longitude of planet c's periastron is slower than that of planet d's periastron except for 
the occasional regressions. Embedded within such a disk, a pair of planets with the present 
observed values of x and rj would not follow a closed hbrating track. When the nebula edge 
is at 5AU (panel 2), the two eigenfrequencies of the system become closer to each other. The 
eccentricities of planets modulate with large amplitudes. The magnitude of x librates on a 
closed equi-Hamiltonian track, about a central point at x ~ 3.3 and r] — 0. 

From the estimation of the equation (D2), we find that Ai — 1 + AN — 7^ vanishes 
when the disk is retreated to Tedge = 5.5AU, namely, secular resonance occurs at Tedge = 
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5.5AU. The panel (3) represents the epoch of secular resonance passage. At this stage, the 
precession frequencies of the longitudes of periastrons nearly match with each other {i.e., 
?7 ~ compared with x). Consequently, rj varies slowly (with finite values other than or tt) 
in comparison with the non negligible evolution of x such that the eccentricities can change 
significantly. The equi-Hamiltonian map is locally symmetric about the lines of 77 = — tt, 
— 7r/2, 0, 7i/2, and vr and that of a; = 7^ ~2. 

When the nebula edge is beyond 5.5AU, the precession speed of longitude of planet c's 
periastron is faster than that of planet d's pcriastron (except for the occasional regressions). 
Panel (4) shows the case that the nebula edge has retreated to 7AU. The centers of the 
closed libration tracks are reversed from the situation in panels (1) and (2). The planetary 
configuration with Cc > and Wc ~ tu^ in panel (2) becomes that with < and zUc ~ zu^ 
in panel (4). When the entire protoplanetary disk is depleted (panel (5)), stable orbits hbrate 
around {x — 0.5, 77 = 0) or (a; = 7, = ±7r). The shape of equi-Hamiltonian map is turned 
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upside down with x — . Figure 3 is the time evolution of the equi-Hamiltonian map for 
the uniform-depletion prescription with rcdge = 4.5AU. In this model, the secular resonance 
occurs when So(t)/So(t = 0) = 0.5. Although the time of the secular-resonance passage 
is slightly modified from that obtained with the model of inside-out depletion, the general 
evolutionary pattern is not qualitatively changed. 

We also numerically integrate the orbits of planets c and d using the the same inside-out 
disk depletion prescription as above {i.e.. we arbitrarily set the initial disk edge to be at 4.2 
AU and specify its retreating speed to be lO^^AU year^^). We consider two sets of initial 
conditions. Figures 4a and c illustrate the evolution of a model with initial values of a; = 8 
and T] = 7r/4. We illustrate the orbital evolution of planets c and d with the {x-rj) diagram 
(see Fig. 4a) and with Cc and as a function of time (see Fig. 4c). These results show that 
the planets' orbits initially librate relative to each other. The planets remain on librating 
closed tracks in the {x-r]) diagram after the disk is totally depleted. But, as a consequence 
of the secular-resonance passage, the eccentricity of planet c becomes smaller than that of 
the planet d. The time-averaged eccentricities after the disk depletion are (cc) — 0.12 and 
{ea) ~ 0.2. 

We also consider a second set of initial conditions with x = 8 {Cc = 0.4, = 0.05) and 
77 = 37r/2. In this case, the planets' orbits initially circulate relative to each other, but, they 
become trapped in librating orbits during the passage through the secular resonance (see 
Fig. 4b). The mean eccentricities after the total depletion of the disk are very similar to 
that obtained in the case of panel (a), because they are determined by the conservation of 
angular momentum (eq. [D6]). 

We also carried out numerical integration of the full equations with the uniform-disk- 
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depletion prescription with t^N — 10^ years, redge = 4.5AU. When the same initial condition 
is used {i.e. x = 8 and rj = 7r/4) were used, the pattern of the orbital evolution with this 
prescription (see Fig. 5a) is not significantly changed from that in Figure 4a. In panels 
b and c, we show the evolution from (a; = 8,?7 = vr) and (x = 6.6,77 = 0). The pair of 
planets in panel b has an initially wide open circulating track. During the disk depletion, 
they temporarily enter into closed librating track. But, their orbital configuration becomes 
wide open again when most of the disk material is depleted. The results of these numerical 
calculations show that as long as the depletion time scale of the disk is longer than the 
oscillation period of eccentricity, the planets' orbits evolve adiabatically. Those systems 
with librating orbits in (x-rj) diagram prior to the disk depletion usually remain on the closed 
librating tracks after disk is totally depleted despite large changes in magnitude of x as in the 
case of Ups And (see below). Those pairs of planets with the widely open circulating tracks 
initially generally remain on open circulating tracks after the disk depletion is completed 
as in the case of HD 168443 (see below). Under some circumstances, it is also possible for 
planets with marginally circulating/librating initial orbits to undergo transition after the 
disk depletion. 



4.2. Orbital Evolution of Circulating Planets Around HD168443 

Next we consider the case of the planetary system around HD 168443 for which we 
assume the mass ratio to be Mb/Mc = 0.451. Figure 6 shows the time evolution of the equi- 
Hamiltonian contour map of the HD 168443 system. The nebula edge retreats from inside 
to outside (from panel 1 to panel 5). The secular resonance occurs when AA^ = —0.86 with 
Tedge = 7.2AU. Similar to the case of the planetary system around Upsilon Andromeda, the 
librating orbit with I77I < 7r/2 are confined in regions with > Cc prior to disk depletion. 
But, in the case of HD168443 system, the domain of closed librating tracks is small. (In 
contrast, for small Mh/Mc{ab/acY^'^, the domain of open circulating tracks is large as is 
the case for planets c and d around Ups And). The domain of closed librating tracks with 
\ri\ < 71/ 2 moves downward in the Hamiltonian map as disk depletion proceeds. If planets 
of HD 168443 initially follow open circulating tracks prior to the disk depletion, they would 
remain on the open circulating tracks after the disk depletion. In this case, during the onset 
of secular resonance, planets b and c may briefly attain a librating track before moving onto 
an open circulating track as the disk continues to deplete. Throughout the epoch of disk 
depletion, the planets' longitudes of periastrons are widely separated such that httle angular 
momentum may be exchanged between them. Consequently, the net change of eccentricities 
is limited. 
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Figures 7 show the numerically computed orbital evolution. The disk potential is calcu- 
lated with the uniform depletion prescription in which we specify ^aat = 10^ years. Figures 
7a and b represent the first model with x = 6 and 77 = initially and Figures 7c and d 
show the second model with x = 1 and rj = initially. For the first model, the eccentricity 
of outer planet c becomes larger than that of planet b after the disk depletion in contrast 
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to the observed eccentricities of planets b and c. With the initial x smaller than 7^ , the 
results for the second model (in panel c) reproduce the present eccentricities of both plan- 
ets. (For presentation simplicity, the notation of 7^ = {ai,/ a^^^'^{Mi,/M^ which represents 
planets b and c.) The extent of the open circulating track is large both before and after the 
disk depletion. During the passage through the secular resonance, the planets' orbits are 
temporarily trapped in the libration cycles around 77 = tt. But the two planets break free 
from their librational state to follow a circulation path after the disk is totally depleted. 



4.3. Dependence on Disk Model 

4..3.I. Single External Disks 

The initial eccentricities of planets, prior to the disk depletion, can be inferred from the 
conservation of the total angular momentum and energy of the planets. We represent the 
mean eccentricities (< Cc > and < >) of the planets c and d around Ups And by their 
values near the libration center {i.e. with x — x^, and = 0). In Figure 8a, we plot the 
values of < Cc > and < > for various Eq and Tedge- The eccentricities < ec > and < > 
are shown by filled circle for six different pre-depletion disk surface densities (5sini, 4sini, 
3sini, 2sini, Isini, and zero times that of the minimum mass model at Tedge where zero 
corresponds to the present state). The model with k = 3/2 has the same power-law index for 
the S distribution as that of the minimum mass nebula model so that these scaling factors 
for S(redge) also represent that for the disk mass. 

We also adopt three different pre-depletion locations of the edge (log(redge/ad) = 0.14, 
0.2, 0.3). The pre-depletion locations of the edge correspond to 3.53, 4.06, and 5.11 AU, 
respectively. The apogee distance of planet d (0^(1 + e^)) is unlikely to exceed the inner 
edge of the outer disk (^cdgc) because its eccentricity may be effectively damped and the disk 
may be strongly perturbed (Papaloizou et al. 2001). From this requirement, we note that 
fedge > 3.5AU for the present value of = 0.35. Prior to the disk depletion, it is possible 
that e'^ may be close to zero. Nevertheless, the tidal torque of planet d can induce the 
formation of a gap with a width which is Ar ~ -\/T2i?Roche where i^Roche = {Md/SM^Y^^ad — 
O.S5ad is its Roche radius (Bryden et al. 1999). Thus, we adopt the minimum value of 
'^edge Od -I- Ar ~ 3.5AU (see top panel of Fig. 8a). 
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The solid line shows an ellipse which is dehneared by the conservation of total angular 
momentum (see eq. [D6] in Appendix). To second order in eccentricity, the conservation of 
the total angular momentum (eq. [D4] in Appendix) implies that 



where and e'^ are the pre-depletion eccentricities of planets c and d respectively. The value 
of 7d ~ 0.286 in the case of planets around Ups And. A maximum value of = 0.37 is 
attained for a minimum value of = 0. But, our results in Figure 4 indicate that x decreases 
during the disk depletion so that > Cc and the maximum value of is unattainable. 

When the disk edge is relatively close to the planet (top panel), > 0.65 and e'^ < 0.1 
if the pre-depletion disk mass is greater than 2 sin i times that of the minimum mass nebula 
model. But for larger values of redge, the pre-depletion eccentricities approach their present 
values as the contribution of the disk mass to the total gravity becomes vanishingly small. 

The secular resonance occurs at ~ 2 in the case of Ups And. As long as there 
was the disk which corresponds to a pre-depletion Sc > 0.51, the planets pass through the 
secular resonance. Note that with an arbitrarily large disk mass, vanishes while attains 
a maximum value of 0.69. Even under these extreme condition, the orbits of the two planets 
do not cross and remain stable. During the disk depletion, the eccentricities of planet c and 
d coincide with each other at about 0.33. For the pre-depletion eccentricity of planet c to 
be > 0.33 (which is generally the case with the exception of the three least massive initial 
disk models) , the inversion of eccentricity ratio between the planet c and d around Ups And 
occurs. 

Because the influence of the disk on the planet is mainly due to the gas near the rdisk, 
the above results are insensitive to changes to the surface density distribution at large radii. 
In Figure 8b, we set redge = 4.057AU (log(redge/ad) = 0.2) and considered three sets of S 
distributions with k — 2, k — 1, and k — 1/2 (Fig. 8b). We adopt six sets of values of the 
surface density S (redge) of each model at Tedge- For disk with similar surface density scaling 
factor, there are no significant differences between three panels. The detailed disk structure 
is not important for the evolution of the planetary eccentricities in our model. 
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4-3.2. Interplanetary Rings and Outer Disks 



After the emergence of multiple planets, gap formation effectively cuts off gas flow from 
the external disk across the outermost planet's orbit. Since both dynamical and viscous 
time scales are rapidly increasing function of the disk radius, the disk interior to the gap 
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is depleted rapidly (Lin & Papaloizou 1986b). Due to a shepherding effect induced by the 
multiple planets (Goldreich & Tremaine 1979), the gaseous interplanetary rings may be 
preserved for a much longer time scale than the internal disk. Nevertheless, the dissipation 
of the competing tidal perturbations by the planets can lead to ring dispersal well before the 
depletion of the external disk (Bryden et al. 1999). 

In the present context, the evolution of the eccentricity ratio (xmo) "the libration center 
around rj = depends intricately on the depletion pattern of both the ring and disk as well 
as the locations of the edges of the ring and disk. For example, the value of x^o is smaller 
for interplanetary rings with inner edge closer to the inner planet. In contrast, the value 
of Xmo is larger for interplanetary rings with outer edge dout closer to the outer planet. The 
balance between these effects determines the evolution of the eccentricity ratio. In order to 
illustrate these dependence, we consider the epoch when planets c and d around Ups And 
were surrounded by an interplanetary ring and an external disk. We adopt a surface density 
distribution which is 5 x sini that of the minimum mass solar model. A hole is cut out 
interior to rfin. An annular strip is also removed between dont and redge- By fixing the inner 
edge of the external disk to be redge = 4 AU, we compute the value of function of 

din (see Fig. 9). In these calculations, we used eight values of the outer edge of the ring dout, 
ranging from 1.2 to 1.9 AU in intervals of 0.1 AU. 

The results in Figure 9 indicate that when rfout < 1.66AU, the gravitational perturbation 
of the ring on the inner planet always surpasses that on the outer planet. Consequently, the 
value of Xmo is reduced below that (ajdisk ~ 12) in the cxtcrnal-disk-only {i.e. ringless) models. 
In these small dout approaches to 12 for relatively large din as the ring's influence 

weakens. But, in the limit dout >1.66AU, Xmo becomes larger than x^isk for relatively large 
din- For example, when d^nt =1.8AU, x^o becomes larger than Xdisk at din >1.5AU because 
the ring's perturbation on the outer planet is enhanced. Also in this large dout limit, 
becomes smaller than current observed eccentricity ratio of 0.66 at din <1.33AU. For all 
values of Oout, the pre-depletion values of Xmo in this model is always larger than the current 
observed eccentricity ratio. 

In Figure 10, we show the evolution of equi-Hamiltonian contours in a model with both 
ring between planets c and d and an external disk beyond planet d. We adopt a prescription 
in which the ring is depleted uniformly prior to the uniform depletion of the external disk. 
For our model, we specify \og{din/ ac) = log(ad/(iout)=0.2 for the ring and log(redge/ad) =0.2 
for the external disk. Our numerical results indicate that Xmo is initially smaller than current 
observed eccentricity ratio (panel 1). The ring proceeds to be depleted during the epochs 
represented by panels (2) and (3). When the ring mass is reduced to half of its original value, 
Xmo become ~ 2 (panel (2)), i.e. the system passes through a state of secular resonance. The 
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entire ring is depleted while the external disk is preserved during the epoch represented by 
panel (3). After that, the contours again evolve through a state of secular resonance (when 
1/3 of the external disk is depleted in panel 4) to current values (in panel 5) as in the case 
of Figure 3. 

Although the same pattern of evolution occurs in the case of planets around HD168443, 

there is little change in the eccentricities of either planets since the contour of their eccen- 
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tricity ratio always retains the values around 7^ . We realize that the evolution of the 
interplanetary ring is determined by the process of gap formation and planet-disk tidal in- 
teraction, its influence on the eccentricity evolution of the planets should be analyzed along 
with its non axisymmetric contributions. The results of such analyses will be presented 
elsewhere. 



5. IMPLICATIONS ON THE PROTOPLANETS' ORBITAL 
ECCENTRICITY PRIOR TO THE DEPLETION OF THEIR NASCENT 

DISKS 

5.1. Planets around Ups And 

First, we consider the evolution of planets c and d around Ups And. We infer from their 
present librating orbits and the results in FigTue 2 that these planets followed a close librating 
track prior to the disk depletion. Our results also indicate that when the disk dominated the 
planets' apsidal precession, the domain of close librating tracks with \ri\ < 7r/2 resided in the 
parameter region with Cc > e^. Thus, the present observed librating orbits of planets c and 
d would be attainable if Cc > prior to the disk depletion. Since Cc < today, the passage 
through the secular resonance must have led to the inversion of this ratio. In Figure 11, we 
map out, in the (x-r)) plane, the most likely initial domain which may have led to the present 
planetary orbits in the Ups And system (panel a). The periapses of those pairs of planets 
which initially occupied region 1 maintain their near alignment despite the large changes in 
X during the disk depletion. These systems enter the domain of closed librating tracks in 
the [x-f]) plane after the disk is completely depleted. Their orbits become very similar to 
the present orbits of planets c and d. Planets in region 2 have open circulating tracks prior 
to the disk depletion. Their orbits remain on circulating tracks after the disk depletion. 
Among those pairs of planets which originated from both regions 1 and 2 {x> ^/^^ ~ 2), 
the eccentricity ratio is reduced after the disk depletion. Thus, we infer those pairs of planets 
which currently occupy the region x{— ec/e^) < 1 are originated from the domain x > 2.5 
prior to the disk depletion. Those pairs of planets with x < 2 initially cannot attain the 
present observed values of x(~ 0.6) because the magnitude of their x increases during the 
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disk depletion. 

In the absence of a strong non axisymmetric torque associated with the planet-disk 
interaction, secular interaction between the planets does not alter the total angular momen- 
tum and energy of the system. In contrast to the eccentricity ratio (x) between planets c 
and d, the magnitude of eccentricities cannot be excited by the sweeping secular resonances 
alone. The equation (17) indicates that if both planets had nearly circular orbits initially, 
their eccentricity would remain small after the disk is completely depleted. Thus, around 
Ups And, at least planet c must have acquired some initial eccentricity prior to the epoch of 
planet formation. 

Up to now, we have neglected the effect of non axisymmetric torque associated with the 
planet-disk interaction. During their formation, the eccentricities of isolated protoplanets are 
damped and excited as a result of their interaction with their nascent disks at the corotation 
and Lindblad resonances respectively. For protoplanets with masses comparable to that of 
Jupiter, the damping effect of the corotation resonances is stronger than the excitation effect 
of the Lindblad resonances (Goldreich & Tremaine 1980). But protoplanets with masses an 
order of magnitude larger than that of Jupiter may open relatively wide gaps in protostellar 
disks. In this hmit, the protoplanets' corotation resonances may be cleared of disk gas such 
that their eccentricity may be excited (Artymowicz 1993; Papaloizou, et al. 2001). For 
modest inclination angle between the their orbital plane and the line of sight, neither planet 
c nor d have sufficiently large mass to clear a wide gap. But, planet c is likely to form prior to 
planet d. The tidal torque of planet d may truncate the disk well beyond the region which 
contains any corotation resonance for planet c so that the latter's eccentricity is excited 
whereas that of planet d is damped. 

The actual mass of planets c and d depends on the value of sin i. If i < 20°, the mass 
of planet d would be sufficiently large for its interaction with the disk to excite eccentricity. 
The constraint on the inclination angle inferred from the dynamical stability of the planetary 
system around Ups And indicate that i > 12 — 45° (Lissauer & Rivera 2001; Ito & Miyama 
2001). In a multiple system, however, the rate of angular momentum transfer between planet 
d and the outer disk may be faster than that between it and planet c. We shall consider in 
a future contribution whether planet c will be able to attain a larger initial eccentricity. 

It is also possible that the planets around Ups And was closer to each other in the past 
and dynamical instability may have led them to undergo orbit crossing (e.g.. Chambers, 
Wetherill, & Boss 1996). Rasio et al. (1996) suggested that such a process may have 
led to the scattering of a planet close to the surface of its host star and subsequent tidal 
circularization may have caused it to become a short period planet. Around Ups And, 
regardless whether such a process can lead to the formation and survival of the short-period 
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planet b, planetary scattering may also have excited the eccentricity and modified the semi- 
major axes of both planets c and d ( Weidenschilling & Marzari 1996, Lin & Ida 1997). 
During such encounters, planets with lower mass generally acquires larger eccentricity as a 
result of the conservation of the system's total energy and angular momentum. Although 
these speculations provide arguments for planet c to have initially acquired more eccentricity 
than planet d, rigorous analysis is needed to quantitative verify these conjectures. 



5.2. Planets around HD168443 

In the case of planets around HD168443, the present orbit follows an open circulating 
track. In Figure lib, we outline the necessary initial domain in the {x-r]) diagram, which 
would lead planets b and c to acquire their present orbital parameters. Planets which 
were initially in the region 2, would enter into open circulating tracks after their nascent 
disk is completely depleted. But their eccentricity ratio would not generally match those 
inferred from the observations of HD168443. The most likely initial domain of planets around 
HD 168443 (with their current dynamical properties) is mapped out as region 1. During the 
passage through secular resonance, the eccentricity of planet c decreases while that of planet 
b increases slightly. The geometry of [x-r]) diagram is reversed during the passage through 
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the secular resonances, i.e., at the point x = . This small change is due to the present 
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value of X being close to that of 7^ in the HD168443 system. In contrast, present x of 
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planets around Upsilon Andromeda is much smaller than 7^ 



5.3. Other Multiple-planet Systems 

We comment here about another system of planets around HD74156 (Table 1(c)). In 
this system, the outer planet is more massive than the inner planet. Since their orbits 
are circulating (similar to those of the planets around HD168443), the evolution of this 
system during the disk depletion would be very similar to that described in §4.2 for the 
planets around HD168443. However, the present value of x is not close to that of 7^^^^. The 
inversion of the eccentricity ratio is expected to occur in the case of planets around HD74156 
(similar to those of the planets around Ups And). 

Our own solar system also contains multiple gaseous giant planets. In this case, the 
inner planet, Jupiter, is more massive than the outer planet, Saturn. Their longitudes of 
periastron are not aligned with each other. The ratio of the Jovian eccentricity to the 
Saturnian eccentricity (ej/es) varies between the range of 0.3 and 5 with time. In this 
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case, the orbits of Jupiter and Saturn do not pass through the secular resonance during the 
depletion of the primordial solar nebula. Therefore, the initial orbits attained by Jupiter 
and Saturn during their formation would not be significantly changed as a consequence of 
the disk depletion. This evolution pattern is similar to that of the two-planet system around 
47 UMa. Those two planets have masses 2.0 and 0.67 Mj and semi-major axes 2.1 AU and 
3.8AU, respectively (Table 1(d)). 

In general, those systems where the outer planets are more massive than the inner 
planets (such as the planetary systems around Ups And, HD168443, and HD74156), the 
depletion of the outer disks induces a secular resonance which sweeps through both planets. 
But in the opposite situation where the inner planets are more massive (more precisely, where 
7^ ' < 1) such as the planetary systems around the Sun and 47 UMa, the depletion of the 
disk does not lead to the secular resonances between the planets. Thus, the changes of the 
eccentricities are not so large. The planets around HD12661 (Table 1(e)) have librating orbits 
with 77 ~ TT. During the disk depletion, the secular resonance passes through this system 
{i.e. Xm attains a value of 7^^^^) because its 7^^^^ > 1 (the inner planet is massive than 
the outer planet). Today, the value of Xm around the libration centers at tt is slightly larger 
than 7(7^^^- Thus, the pre-depletion value of x^ must be smaller than 7^^*^^. Nevertheless, 
the magnitude of the eccentricity change is relatively small because the present value of x^ 
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is close to that of 7^ 



6. SUMMARY AND DISCUSSIONS 

In this paper, we consider the secular interaction between pairs of coplanar planets. 
Provided these planets are not locked in a mean motion resonance, their secular interaction 
leads to angular momentum exchange without energy transfer between them. Consequently, 
both planets undergo apsidal precession and eccentricity modulation about some mean val- 
ues. Thus, dynamical evolution of these systems of planets is best illustrated by the ratio of 
their eccentricity x and the relative longitude of their periastrons 77. For example, based on 
their observed orbital elements, we find that gravitational perturbation between the outer 
two planets around Ups And causes both 7] and x to modulate, with small amplitude, along 
a close librating track around 77 = and x ~ 0.66 (see Fig. 1). 

Although the mean eccentricities of these planets are constant today, they may have 
undergone significant changes during the epoch when their nascent disks were depleted. In 
order to extrapolate the initial dynamical properties of planetary systems, we consider the 
precession caused by the self gravity of these disks. We show that contribution from the 
disks initially dominates the precession frequency of the planets and regulates the angular 
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momentum exchange rate between them. As their influence fades during the disks' depletion, 
the precession eigenfrequencies of the planets' longitude of periastron pass through a state 
of secular resonance. 

As a result of this transition, the magnitude of both planets' eccentricities changes 
with their ratio inverted. But, provided the disk depletion proceeds on a time scale long 
compared with the precession frequencies associated with the secular interaction between the 
planets, the libration amplitude of t] is not significantly changed due to the preservation of an 
adiabatic invariance. Thus, if a pair of planets started out with a close librating orbit prior to 
the disk depletion, they would generally retain the librating nature of their secular interaction 
after the disk depletion. Planets with open circulating tracks initially would generally evolve 
to attain other open circulating tracks after the disk depletion. However, the circulating 
track may temporarily enter into closed librating track near the secular resonance. Planets 
with marginally circulating orbits (close to librating orbits) prior to the disk depletion may 
sometimes attain librating orbits after the disk depletion provided the difference between the 
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libration and resonant centers {i.e. \xm — Id I) reduced during the transition. Inversely, 
planets with marginally librating orbits prior to the disk depletion may also sometimes attain 
circulating orbits after the disk depletion provided the difference between the libration and 
resonant centers is increased. 

Using these results we infer the kinematic properties of planetary systems around both 
Ups And and HD168443. During the passage of the secular resonances, the libration time 
scales of the planets around these two stars are 2 x 10^ and 10"'' years, respectively. Provided 
the external axisymmetric disk is depleted on a time scale longer than these libration time 
scales, transition through the secular resonance is adiabatic. 

In this contribution, we show how the process of disk depletion may lead to signifi- 
cant modifications in the eccentricity distribution within some extrasolar multiple planetary 
systems. These changes are the result of secular perturbation the planets exerted on each 
other. But planetary secular perturbation and the axisymmetric potential of evolving disks 
alone do not lead to eccentricity excitation. We have not quantitatively addressed the issue 
of initial eccentricity excitation among isolated planets or at least one planet in a multiple 
planetary system. We will discuss these problems in the future. 
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A. PRECESSION DUE TO THE DISKS' GRAVITY 

In this section, we evaluate AA^ in equation (9). After protoplanets have acquired a 
sufficiently large mass, they induce the formation of a gap near the orbit through tidal 
interaction with the disk (Goldreich & Tremaine 1980; Lin & Papaloizou 1980, 1993). While 
the disk material can no longer be rapidly accreted onto the protoplanet (Bryden et al. 1999), 
it continues to contribute to the gravitational potential. Although the disk does not have 
a sufficiently large mass to modify the angular frequency of the planet's orbit, its gravity 
can induce precession (Ward 1981). The disk's contribution to the gravitational force in the 
radial direction at r is 

where S is the surface density of the disk, ^ is the potential, y = r'/r, ip = 2|/^/^/(l + y), 
K{ijj) and £'('^) are elliptic integrals of first and second kind, and G is the gravitational 
constant (eq. 2-146 in Binney & Tremaine 1987). The associated precession frequencies for 
planets c and d can be obtained from 

N,, = (Qe,. - = ^ (^Ir'F.) ^ ^ , (A2) 

where Qc,d and k^j are respectively the angular and epicycle frequencies at semi-major axis 
of planets c and d (a^j). Depending on the S distribution, Nc^^ ~ 0{nf.,dtcM£,/ M^) where 
Md = Ea^ ,^ is the characteristic disk mass at r = ac,d respectively. From equation (A2), we 
find 

with its sign determined by the S distribution. For disks with finite S only outside the 
planets' orbit, the addition of the disk self gravity leads to both Nc and Nd to be positive. 
Since planet d is closer to the external disk, N^, > so that AA^ is negative. For disks 
which engulf both planets' orbits, AA can still be negative provide E does not decrease too 
rapidly with r. For example, if we use the minimum mass nebula model (Hayashi et al. 
1985) as a fiducial prescription for E such that 

S = 1.7 X 10^(r/lAU)"^-Vm"^ (A4) 



AA is negative. 



B. ORBITS OF LIBRATING AND CIRCULATING PLANETS 



Using linearized equations (8) and (9), we consider here three famihes of orbits in systems 
which contain two planets. We assume that the two planets orbit on the same plane as the 
disk. 



As an example, we consider the planets c and d around Ups And. For planets c and d 
around Ups And, C = —0.4, 7^ = 0.29. The longitudes of these two planets are aligned such 
that their 77 ~ today. The present values of a; = 0.66 and tc = l-l x 10^ years. With these 
values, we show below that secular interaction between the two planets induces x and rj to 
undergo small amplitude oscillation around stationary points. These points correspond to 
extrema for x and 77. Prom the linearized equations (8) and (9) of secular motion, we find 
that these extrema occurs at 77 = and tt. 

For 1] = 0, a stationary solution can be obtained for x — Xm{> 0) where x^ is a root of 
the quadratic equation 



With X = Xm and rj = 0, longitudes of the periapse of both planets are aligned and precess 

together, even though gi ^ g^- Consequently, there is no angular momentum exchange and 
the eccentricities of both planets are conserved. For solutions near the stationary point with 
7] = TT, Xm can be obtained from a similar equation where 



In the {x-Tj) diagram, these orbits follow the equi Hamiltonian contours (see Appendix 
D). Angular momentum is being continually exchanged between the planets, leading to an 
eccentricity modulation. However, the libration is centered around 77 = or ±7r where there 
is no net torque between the two planets. Although the range of 77 modulation (A77) is 
limited, the two planets are not in the center of secular resonance because 77 modulates with 
a frequency uj which does not vanish with arbitrarily small A77. 

We now consider a small perturbation about the stationary point. For small values of 
77 and e = x/ Xm — ^1 equations (8) and (9) can be linearized in 77 and e such that 



B.l. The Orbits of Librating Planets 






(B2) 



de 



X, 



C 



(l + 7dX^)77, 



(B3) 
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Combining these two equations, we find the solution 

e = eosin(a;T), 77 = 77oCOs(a;T), with 770 = ±eo, (B5) 
where the osciUation frequency is 

^ = — (1 + l<i^l) , (B6) 

where the signs of the expression for rjQ and u are the same. Near the stationary point 
X = Xm and 77 = 0, a; is also approximately the difference of the precession frequencies of the 
two planets because zu for both precess in the same direction. 



B.2. Orbit of Circulating Planets 

As another example, we use the planets around HD 168443. For planets b and c around 
HD168443, C = —0.13, 7^ = 0.14, and tc = 1.9 x 10^ years. For notational consistency, we 
use the subscription for the inner planet is c and that for the outer planet is d. With these 
parameters, AN — for Xm = 0.15 (or 46) and = (or tt). The observed value of a; = 2.65 
and 77 = 1.92 which is far from the stationary point x — Xm and 77 = (or tt). With these 
parameters, \Ai\ ^ IA2I (see eqs. [7] and [9]), we expect the orbits of planets b and c to 
circulate with respect to each other. In this case, the hnearized equations (B3) and (B4) are 
no longer applicable. The small magnitude of C in equation (8) implies that the modulation 
amplitude of x is small compared with its average value. The dominance of the term of 
1 — + AA^ in equation (9) implies that, to zeroth order in r, a; ~ Xq and 77 ~ (1 — 7^)7" 
where Xq is a constant. Expansion to the next order in r yields the solution that 

77 (1 - 7. + mr + - ,/^~TpAn ^^"(^ - + ^^)^' (B^) 

Xq (1 - 7d + AiV) 

X ~ xo - ^— — cos(l - 7d + AiV)r, (B8) 

1 - 7d + A7V 

so that the circulation period is ti ~ 27ric/(l — 7d + AN) ~ 1.4 x 10"^ years for AN — 0. 



B.3. Non Periodic Orbits Near Secular Resonance 

In addition to the librating and circulating orbits, there is another family of non periodic 
orbits for the case when \Ai\ < \A2\. In the limit that both x and 77 have values which are 
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very different from Xm and 77=0 (or tt) , r] would evolve rapidly to a phase angle 



such that drj/dr is, to within the first order approximation, reduced to zero (see eq. [9]) and 
the two planets become phase locked. Within the same order approximation, x becomes Xi 
and the evolution of x is determined by 

provided \Ai\ remains smaller than \A2\. (In Appendix C, we consider the depletion of the 
disk during which the magnitude of Ai evolves with time.) 

The plus and minus signs on the right hand side of equation (BIO) correspond to < 
rji < 71 and n < rji < 271 respectively. Increases/decreases in x are equivalent to eccentricity 
excitation/ damping for body c. As a consequence of x evolution, rj also adjust to ~ 771 such 
that the evolution of the planets' orbit is non periodic. Nevertheless, during the raise/decline 
in X, the two planets may also undergo small amplitude {x2, 772) libration around the evolving 
guiding center at xi and 771. The evolution of X2 and 772 can be determined from the next 
order expansion of equations (8) and (9) such that 

^^B,7J2 + B2X2, (Bll) 

^ = Bo{t - To) + BsX2 + B^ri2. (B12) 
ar 

The coefficients are 

Bo = AN, 5i = C(l + 7dX?)cos77i, 

B2 = 2C7da;isin77i, B3 = -{C/xl){l + 'ydX^cosrii, 

- -{C/xi)il-^dxi)smr]i, (B13) 

where AA^ = dAN/dr and the fiducial dimensionless time tq corresponds to the instance 
when the two planets first enter into a phase lock, i.e. when their x — Xi and r] — rji. The 
solution of equations (Bll) and (B12) can be expressed as 

X2 = X2oexp{a;2o(T-To)} + X2iexp{a;2i(T-To)} + -B5T + S6, 

772 = ^ [((^20 - -B2)a;2oexpa;2o(r - To) + (ci;2i - -B2)a;2iexpa;2i(T - To)] 

+^ [^5(1 - B2T) - B2Be] , (B14) 
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where X20 and X21 are the amphtudes of oscillations and 



B2B4 — B1B3 



BqBi 



and Bq = 



{B2 + Bi)B^ — BqBiTq 
B2B4 — B1B3 



(B15) 



with the growth rate 




equations (Bll) and (B12) becomes slightly modified from that in equations (B14). Also 
note that for r/i = or tt, so that equations (Bll) and (B12) reduce to equations (B3) 
and (B4) respectively. Also, UJ2 becomes purely imaginary {±iu) so all the orbits undergo 
libration about the stationary point. But, for a more general value of 771, UJ2 is complex, 
and oscillations about the evolving xi values either grow or decay depending on the value 
of rji. At the center of resonance where Ai = and rji — 7r/2 so that 002 — 2(77^X1 or 
{C / Xi) {'ydxl — 1) ■ Both roots are real so that the X2 either converge or diverge (if xi > "y~^^^) 
without any oscillation. Note that oscillation about Xi and rji only occurs if the magnitude 
of UJ2 is larger than B^ such that changes in AA^ and Xi occur in an adiabatic manner. For 
further discussions on the evolution of planetary orbits during the epoch of disk depletion, 
see Appendix C. 



C. THE EVOLUTION OF MULTIPLE-PLANET SYSTEMS DURING THE 



During the depletion of the disk, changes in the magnitude of AA^ can cause Ai to 
change sign. In §4, we presented numerical results of orbital evolution during the epoch 
of disk depletion. Here, we present analytical solutions to interpret the numerical results. 
Following the discussions in Appendix B, we consider three initial orbits separately. 

C.l. The Evolution of Two-Planet Systems with Initially Librating Orbits 

In Appendix B.l, we show that the condition for libration is \Ai\ < \A2\ and eo ^ 1. 
For planets c and d around Ups And, AA^ = today so that Xm ~ —C/ (1 — 7^) ~ 0.5 and 
(J ~ 0.83. Since the observed values of x and rj are close to those of x^ and respectively, 
the small- amplitude oscillator approximation is valid. The necessary condition for libration, 
i.e. \Ai\ < 1^2! is marginally satisfied. From equation (12), we find that the time scale for 
one complete libration cycle today is ti — Tubtc — 27rtc/u; — 8.5 x 10^ years. 



EPOCH OF DISK DEPLETION 
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Prior to its depletion, five times minimum mass nebula model would lead to Mo/Ma — 
Ti{ad)o?d/Md ~ 0.015 and ~ —1.4 (eqs. [A3] and [D2]). In this case, the necessary 
condition for libration is also marginally satisfied and it is possible to find a real positive 
solution to the quadratic equation (Bl). Around the stationary point Xm ~ 6.5, there are also 
solutions which correspond small amplitude libration between planets c and d (|?7| < 7i"/2) 
despite the rapid precession induced by the disk's self gravity. The corresponding uj ~ 0.80 
for such a protoplanetary system is very similar to the present disk- free value of Xm- 

During the depletion of the protostellar disk, the magnitude of AA^ gradually vanishes 
so that both Xm and u changes with time. For —1.5 < AA^ < 0, the condition for libration 
(|y4i| < 1^2!) is satisfied. Provide rub < r^i. Tut < r^^, and < r^^ where 

_ 27r 2ttx _A^ 1 - 7rf + AAT 

UJ C{l + -idX^) Ai AN 



Tx™ = — = - J- Ta,, T^ = - = [ TA„ (02) 



the libration of the planets' orbits is preserved while the librational center and frequency 
evolves adiabatically. 

At the special phase when Ai — 0, the second term in equations (Bl) and (B2) reduces 
to zero and the stationary points with = and rj = tt coincide with each other with 

— 1/2 

Xm = 7d ■ refer to this special configuration as librational degeneracy. In this state, 

— 1/2 

all orbits librate cither around Xm = Id , rj = or n with two separatrix located on the 
r] = n/2 and 37r/2 (see Hamiltonian contours maps). From equation (6), we note that in 

this degenerate state, Agfi2 attains a minimum value with respect to variations in AA". The 

1/2 

oscillation frequency around the stationary points u attains a minimum magnitude 27^' C 
with Xm — 7d (see eq. [12]) which is ~ 1.8 for planets c and d around Ups And. 

Note that in this state of librational degeneracy, gi 7^ g2 in general. However, if planet 
c has negligible mass, Agi2 would vanish and it would enter into a state of exact secular 
resonance and x^ would diverge. This asymptotic state represent the eccentricity excitation 
of asteroids by the Jovian planets. For finite values of 7^, this state of hbrational degeneracy 
can also be considered as a state of secular resonance despite the finite value of Ag'12. The 
secular interaction between the planets forces them to librate so that their relative longitude 
of periapse passage is confined to a limited amplitude. 

We now consider the evolution of the librational amphtude. The action associated with 
the simple oscillator equations (B3) and (B4) 

ede — e^uj / cos^OdO — Tie^uj (C3) 
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is invariant to first order in tya^/ta^ or riib/T(^ (Goldstein 1980), which ever is larger. (For 
values of r^^ and t^, see eqs. [CI] and [C2]). Based on this invariance, we can also determine 
the changes in eo which is oc For planets c and d around Ups And, we adopted a 

model in which uo ~ 0.8 both during the embedded phases and after the disk is completely 
depleted. Thus, we expect both eo and r^o to be remain invariant despite large changes in 
Xm- The numerical results in §4.1 is in agreement with this analytic derivation. 

But, the condition for adiabatic evolution may not always be satisfied especially when 

— 1/2 

the two planets pass through a phase of secular resonance with Xm = 7^ when the libration 
period takes the maximum value ti = -T^tcl {C%1'^) = 1.7 x 10^ years and lu attains a 
minimum value of ±26*7^/^ = ^0.43 (compared with its initial value of ~ 0.8). During this 
phase, the libration time scale Tub reaches a maximum. Although at the center of secular 
resonance, ta^ vanishes with Ai, the more appropriate value for it is ~ AN/ AN which 
remains finite. Nevertheless, it is possible for rub to be larger than ta^. Note that 
reduces to ~ ±2/ AN in the resonance. In accordance with equation (C2), t^^ diverges. But, 
the characteristic time scale for u! evolution based on a second order expansion remains finite. 
The numerical results in §4.1 show a slight increase in eo with a ~ 6% reduction in rjQ. In 
contrast, a 40 % increase in eo is inferred from equation (C3). Nevertheless, for planets c 
and d around Ups And, the duration of non adiabatic evolution is brief so that the small 
amplitude libration around x — Xm and rj = is preserved despite the large net changes in 
the magnitude of Xm- 

Note that wc have only consider the axisymmetric contribution of the disk to the total 
potential. Therefore the total angular momentum of the two planets is conserved during 
the disk depletion even though there is a net transfer of angular momentum between them 
which leads to the decline in Cc and the growth in e^. 

C.2. The Evolution of Two-Planet Systems with Initially Circulating Orbits 1: 

Asteroids 

We now consider the orbital evolution of multiple- planet with initial values of |y4i| > 
\A2\. In Appendix B.2, we showed that these systems have initially circulating orbits. There 
are two possible branches of solutions. First, we consider the limiting case of small 7^ 
(<S x'^) which is appropriate in the context of secular interaction between the asteroids and 
Jupiter for which Mc < M^. 

During the epoch of disk depletion, as the magnitude of Ai in equation (9) is reduced 
below that of ^2, V evolves to rji (see eq. [B9]) which reduces to 771 = cos~^{—Aix/C) in the 
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limit 7d <S x~'^. Since Ai, A2, and C are negative, rji — ir when the planets first enter into 
a secular resonance. Also, equation (BIO) implies that Xi decreases/increases for rji slightly 
less/greater then tt. 

Let us first consider the case that rji is slightly less than tt. Continual decline in x is 
attainable provided drji/dr < 0. The reduction in x leads to an increase in the magnitude 
of A2 while the disk depletion causes a decrease in that of Ai. Consequently, the magnitude 
of A1/A2 decreases in equation (B9), leading to a reduction in rji. This state of phase lock 
can only be maintained if rj can adjust to 771 and decrease to ~ 7r/2. So long this state of 
resonant phase locking is maintained, the magnitude of x {i.e. Xi) can decrease indefinitely 
When Ai passes through the zero value and becomes positive, x continues to decline. If the 
magnitude of x cannot be reduced as fast as that of AN, \Ai\ > \A2\ and the planets would 
attain circulating orbits. Even in the limit that the magnitude of x can be reduced as fast 
as that of AA^, the planets are unlikely to attain hbrating orbits because their rji ~ 7r/2. 
Eventually, rji terminates its evolution with values slightly less than 7r/2 while the orbits of 
the massless particles are approximately circularized. 

If the planets enter into the secular resonance with rj slightly larger than tt, x would 
increase. Further excitation of x requires drji/dr > 0. In equation (9), the growth in x 

reduces the magnitude of A2 while that of Ai is also decreasing. Phase lock can only be 
maintained if \Ai\ < \A2\ such that x cannot grow faster than C/Ai ~ C/(l — AA^) during 
the epoch of disk depletion. If x grows slower than C/Ai, drj/dr > and rii would increase 
to enhance the x growth rate. The magnitude of rji is determined by the ratio of Ai and A2 
(see eq. [B9]). The evolution time scale of A2 equals to Tx^. Equations (CI) and (C2) imply 
that — in the limit of small 7^^^. Thus, magnitude of A1/A2 is regulated to be less 
than unity as phase lock causes x to increase (see eq. [BIO]). 

As the resonant condition 1 — AA^ ~ is approached, rii increases to 37r/2 and x grows 
at a maximum rate dx/dr = —C. So long this state of resonant phase locking is maintained, 
the magnitude of x can increase indefinitely. On the resonance, gi = §2 such that the growth 
occurs regardless of the value 77. Increases in x is equivalent to eccentricity excitation for 
the less massive body c. When Ai passes through the zero value and becomes positive, x 
can have a sufficiently large value such that Cc may approach and even exceed unity. This 
eccentricity excitation process has been noted by Heppenheimer (1980) as a potential cause 
for the existence of highly eccentric asteroids. 

As X increases in the resonances, there are three eventual outcomes: a) x becomes 
sufficiently large for Cc to exceed unity and body c to escape, b) the planets would attain 
circulating orbits as the magnitude of Ai exceeds that of A2 after |AA^| decreases below 
~ 1 if becomes much smaller than r^j^, and c) in systems with finite 7^, the above 
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approximation become inappropriate when x exceeds 7^ (see below) . 

C.3. The Evolution of Two-Planet Systems with Initially Circulating Orbits 2: 

Two Planets with Comparable Masses 

We now consider the second branch of solutions in which 7^ is finite. These parameters 
are more appropriate for the extrasolar multiple-planet systems in which the masses of the 
planets are comparable. The main difference introduced by the finite value of 7^ is the 

— 1/2 

possibility of a transition as the value of x passes through 7^ . This branch of solution 
is applicable for planets b and c around HD 168443. We use the analytic solutions in this 
subsection to study the numerical results in §4. 

For planets b and c around HD168443, C = -0.13, 7^ = 0.13, and = 1.9 x 10^ years. 

In this system, the value of x,„ depends sensitively on AA^ due to the small magnitude of 
C and 7d. Within the uncertainty of analytic estimate, we estimate M^/Md ~ 0.007 and 
AN ~ —4 for five times minimum mass nebula model. Despite this relatively small value of 
Md{< Mfi < M^,), the corresponding Xm — 180. At this initial stage, \Ai\ > \A2\ and the 
planets follow circulating orbits which is well approximated by equations (B7) and (B8). 

But, during the depletion of the disk, the decline of S reduces the magnitude of AA^ 
which passes through a stage of secular resonance in which AA^ ~ 7^^ — 1 = —0.86 with the 

— 1/2 

corresponding value of — 7^ = 2.6 similar to the present value of x. In this limit, 
I All < 1^421 and the leading-order approximation with the solutions in equations (B7) and 
(B8) are no longer vahd. Since the planets' orbits circulate relative to each other prior to 
entering this stage, the small-oscillation-amplitude approximation which led to the solutions 
in equation (B5) may also be inappropriate. Nevertheless, when the planets first enter into 
the resonance, \Ai\ = IA2I, f]i = tt, and UJ2 in equation (B16) becomes purely imaginary. 
The discussions in §B.3 indicates that during the passage through the resonance, all orbits 
are librating regardless of the amplitude of the libration cycles. 

Within the resonance where \Ai\ < \A2\, the two planets become phase locked such that 
Xi and ?7i evolve in accordance to equations (B9) and (BIO). Depending on the rate of disk 
depletion, there are three potential outcomes for passage through the secular resonance. 1) 
In the limit that r^-^ < rub, the system evolves in an impulsive manner. Such a situation 
is relevant for rapid clearing of disks which contain long-period, low-mass, widely-separated 
planets with large tc- In this case, there is insufficient time for the longitude of periapse to 
significantly circulate before the effect of secular resonance fades with the rapid disk deple- 
tion. During this rapid passage through the secular resonance, both x and rj essentially retain 
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the values they had when they first entered into the resonance {xi and 771). These values of 
Tji and Xi become the initial values of circulating orbits in the new disk-free potential. Since 
rji and Xi are random parameters, the actual variational amplitude of the circulating orbits 
are unpredictable. 

2) For systems such as planets b and c around HD168443, < Tub so that the planets 

adiabatically evolve into a phase lock. The relevant governing equations become equations 
(B9) and (BIO). In contract to the case for the masslcss asteroids, the system first enters or 
finally exits the secular resonances with = if a; > 7^ or tt if a; < 7^ ^^'^ . In the middle of 

—1/2 

the resonant zone where 7^ — 1 — AN = 0, although Xm = Id ' '"^ generally docs not equal 
Xm so that T] = ±7r/2. (Even in the limit that x = Xm at the center of resonances, A1/A2 also 
vanishes at the center of the resonance). Our analysis in §B.3 indicates that, at the center of 
resonance, both roots of 0J2{— 2C^dXi or (C/xi)(l — 7^2;^)) are real so the amplitude of small 
perturbations either converge or diverge (if xi > 7^ ^^^) from the stationary points without 
any oscillation. Note that in contrast to the case of massless asteroids, uj2 changes sign as 

— 1/2 

xi passes through 7^ so the system can neither be damped nor be excited indefinitely. 
This sign reversal of uj2 limits the among of changes in xi. Physically, this limit is due to 
the comparable momentum of inertia carried by planets of similar masses. In such systems, 
the exchange of angular momentum cannot lead to vastly different eccentricity changes for 
the two planets. 

After passing through the center of the resonance, rji continues to evolve towards or 
TT in accordance to equation (B9) while Xi evolves in accordance to equation (BIO). When, 
rji reduces below sin^^O.8, uji once again become complex and libration around xi and rji 
would resume if the magnitude of UJ2 is larger than such that changes in AA^ and xi 
occur in an adiabatic manner. But in equation (B15), the magnitude of is determined by 
A7V or equivalently ta^- Thus, in the modestly rapid depletion case where r^i < the 
difference between xi and x^ at the center of resonance is preserved. For small C magnitude 
(as in the case of HD168442), the width of the libration is small and Xm is also not changed 
significantly during the passage through the resonance. Thus, planets would not be able to 
enter into a librating state and the circulation pattern of their orbits is retained. 

3) Finally, in the slow depletion limit where > <^2^^, the protracted passage through 
the secular resonance can greatly modify both x and r]. As the planets exit the resonance, 
they can become trapped onto librating orbits. Large changes in x implies that at least one 
planet attains a relatively large eccentricity. In relatively compact multiple-planet systems, 
large eccentricity can cause mean motion resonances to overlap and dynamical instability. 
We shall address the issue of dynamical instability elsewhere. 
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D. HAMILTONIAN CONTOURS 

These phase space evolution of the planets' orbits may be compared with equi-Hamiltonian 
contour maps. By expanding the disk potential in equation (14) for planets with eccentric 
orbits and then averaging it over their orbital period, Ward (1981) obtained the disturbing 
function due to the disk potential as follows: 

< .1, . " t A.^^^ (^) . 4,T,. (Dl) 



Using this potential, AA^ of equation (A3) can be written as 



AiV = ^^^.-T.g)"^^ (D2) 



When the disk depletion time scale is longer than characteristic librational time scale as- 
sociated with equations (8) and (9) (time scales are given in Appendix C), AA^ can be 
regarded as quasi static. In this limit, the planets' orbits follow a single trajectory on {x-r]) 
plane depending on the magnitude of their Hamiltonian and angular momentum. Unless the 
orbits of the interacting planets are locked in mean motion resonances, the orbital energy 
and semi-major axes of each planet are conserved to the lowest order. The orbit-averaged 
Hamiltonian (H) can be calculated semi- analytically using the disturbing functions of plan- 
ets (e.g., Brouwer & Clemence 1961) and equation (Dl). Expanding to the second order of 
the eccentricities, we find 

Menial Mdulal UculNc , naalNa , GMcMaUc 



2 2 2itc 'itc 80^ l ' ' J 

(D3) 

From the conservation of the total angular momentum, a relation 

J = McUcal^/l - el MdUdalyJl - = constant (D4) 

holds provided their orbital planes are not inclined from each other. Thus, to second order 
in eccentricity, 

UcalN^ naalNd GMcMdttc / (1) . 2 , 2n oa(2) \ + + 

— + — + ^ {%/2{ec + - 26^/2606^ cos 77| = constant, (D5) 

McUcalel + MdUda^el = constant, (D6) 
for any given values of H and J. Divided equation (D5) with equation (D6), we find that 

x^l + Ae) + 2Cxcosr/ + (7^^ + Nd) 



constant, (D7) 
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and obtain the equi-Hamiltonian contour lines for various values of the constants of integra- 
tion. In equation (D7), the masses of the disk and planets are contained only in the form 
of their ratios. Therefore, as long as M^, ^ Mc^^, the equi-Hamiltonian contours are not af- 
fected by the obliquity i of the system. We show below that other than small deviations due 
to non-linear effects, the planets' orbits closely follow the contours in the equi-Hamiltonian 
map. 
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Fig. 1. — Orbital evolutions of extra solar planets. Points show the orbits obtained from 
numerical simulations. The dotted lines are the equi-Hamiltonian contours of the planetary 
system. The contours are spaced with equal intervals, (a) Orbital evolution of planets c and 
d (points) of Upsilon Andromedae system. Orbits are integrated for 2 x 10^ years starting 
from present orbital elements. The period of the libration is about 6000 years, (b) Orbital 
evolution of planets b and c (points) of HD168443 system. Orbits are integrated for 2.5 x 10^ 
years. The circulation period for rj to evolve from — tt to tt is about 16000 years. 
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Fig. 2. — Equi-Hamiltonian contours of Upsilon Andromeda planetary system in [x-r]) dia- 
gram. The time sequence of nebula depletion proceeds from panel (1) to (5). The nebula 
edge is set to be (1) 4.2AU, (2) 5AU, (3) 5.5AU, (4) 7AU, and (5) infinity (the present 
disk-free planetary system). 
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Fig. 3. — Equi-Hamiltonian contours of Upsilon Andromeda planetary system in [x-r]) dia- 
gram. The time sequence of nebula depletion proceeds from panel (1) to (5). The depletion 
ratio exp(-t/tAiv) is at (1)1, (2) 0.75, (3) 0.5, (4) 0.25, and (5) 0. The disk edge is located 
at 4.5AU. 
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Fig. 4. — Orbital evolutions of planets during the depletion of the protoplanetary disk, (a) 
In this system, planets were initially on a librating track in the [x-r]) diagram and remained 
on a hbrating track with a high after their nascent disk is completely depleted. The initial 
magnitude oi x — 8 and rj — n/A prior to the epoch of the disk depletion, (b) In this system, 
planets were initially on a circulating track in the {x-rj) diagram but entered on a librating 
track after the disk is completely depleted. The initial magnitude of Cc — 0.4, = 0.05, 
and 7] = 37r/2 prior to the disk depletion, (c) The evolution of the eccentricities in the case 
of (a). The solid and dashed lines show the eccentricity of planet c and d, respectively. The 
retreating speed of nebula edge is 10~^(AU year~^), and the nebula edge is 4.2 AU at t = 0. 
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Fig. 5. — The orbital evolutions of planets during the uniform depletion of the protoplanetary 
disk, (a) In this system, planets were initially on a librating orbit. The initial condition is 
the same as the case of Figure 4a. (b) In this system, planets were initially on an open 
circulating track and remained on an open circulating track after the disk is completely 
depleted. The initial magnitude of Cc — 0.4, = 0.05, and r] = n. (c) In this system, the 
coincidence of the pcriastrons is preserved. The magnitude of Cc = 0.45, = 0.08, and 
77 = initially. In these calculations, the nebula edge is set to be 4.5 AU, and t^N = 10^ 
years. 
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Fig. 6. — Equi-Hamiltonian contours of HD 168443 planetary system in the {x-r]) diagram. 
The time sequence of the nebula depletion proceeds from panel (1) to (5). The nebula edge 
is at (1) 5AU, (2) 6.5AU, (3) 7.2AU (secular resonance) (4) 8AU, and (5) infinity (present 
planetary system). 
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Fig. 7. — Orbital evolutions of HD 168443 planets during the uniform depletion of the proto- 
planetary disk. The edge of the disk is located at rcdgc = 5.5AU. The disk depletion timescale 
is t^iv = 10^ years, (a) The orbital evolution in the {x-r]) diagram for the case with the initial 
values of a; = 6 and 77 = 0. The eccentricity ratio for the first t = 4 x 10^ years is denoted 
by filled squares and that for t = 4 x 10^ ~ 7 x 10^ years (around the secular resonance) 
is denoted by open circles and that for i > 7 x 10^ years are shown by small plus symbols, 
(b) The evolution of eccentricities corresponds to the system of planets in panel a. (c) The 
evolution of the eccentricity ratio for the case with the initial values of x = 1 and rj — 0. 
Filled squares, open circles, and small plus show the period of time before 4 x 10^, during 
4 X 10^ ~ 6.2 X 10^, and after 6.2 x 10^ years, respectively. 
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Fig. 8. — Initial eccentricities of planets c and d around Ups And. The eccentricities Cc 
and Cd are shown by filled circle for five deferent initial disk mass. The number along the 
filled circle shows the surface density (e.g., '5' means Ssinix of the minimum mass model). 
(a)The locations of the disk edge are log(rcdgc/ad) = 0.14 (top panel), 0.2 (middle panel), 0.3 
(bottom panel), (b) The surface density gradients are k = 3/2 (top panel), k = 1 (middle 
panel), k — 1/2 (bottom panel). The location of the disk edge is 4.057AU. 
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Fig. 9. — The location of the hbrating center around zero {Xmo) in the case that there is 
a ring between the planets c and d around Ups And. The ring's surface density is set to 
be 5sini times that of the minimum mass model and it is confined between rfin and dout- 
Eight lines for different values of dout are shown. The horizontal line at ~ 12 shows the x^o 
without the disk. The horizontal line at 0.66 shows the current x^o- 
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Fig. 10. — Equi-Hamiltonian contours of Ups And planetary system. The surface densities 
of residual ring and the disk are (1)100% , (2) 53% (ring) and 100% (disk), (3) 0% (ring) and 
100% (disk), (4) 0% (ring) and 35% (disk), and (5)0% of their initial values, respectively. 
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Fig. 11. — (a) A schematic diagram of the [x-r]) plane for the system of planets around Ups 
And. Planets which initially occupy region 1, enter into a librating track in the {x-r]) plane 
after the disk depletion. Their orbits become similar to present orbits of planets c and d. 
For planets which initially occupy region 2, their eccentricity ratio is reduced after the disk 
depletion, (b) A schematic diagram of the [x-r]) plane for the planetary systems around 
HD168443. Planets which initially occupy region 2, enter into the open circulating tracks in 
the [x-f]) plane after the disk depletion. Planets which initially occupy region 1 evolve into 
orbits similar to those observed around HD168443. 
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Table 1. Orbital parameters of extrasolar multiple planets planets 



planets a(AU) e 


ro(rad) 


Msini(Mj) 


T ■ 

pen 


(a) Upsilon Andromedae (M* = I.SMq) 








planet b 0.059 0.01 


5.522 


0.69 


2450000.6383 


pidnet c u.oz( u.zo 






9/1 cin'iQQ n 


pidlltlL LI Zi.UU U.OU 


4 XT A 


A in 


94^1 "lAR Q 


(b) HD168443 (M* = I.OIMq) 








planet b 0.295 0.53 


3.018 


7.73 


2450047.58 


planet c 2.87 0.20 


1.098 


17.15 


2450250.6 


(c) HD74156 (M, = LOSM©) 








planet b 0.276 0.649 


3.206 


1.56 


2451981.4 


planet c 3.47 0.395 


4.189 


>7.5 


2450849 


(d) 47 UMa (M, = IMMq) 








planet b 2.09 0.061 ±0.014 


3.00 ±0.26 


2.54 


2453622±34 


planet c 3.73 O.liO.l 


2.22 ±0.98 


0.76 


2451363.15 ±493 


(e) HD12661 (M, = I.O7M0) 








planet b 0.804 0.35 


5.11 


2.21 


2460208.393 


planet c 2.652 0.11 


2.29 


1.58 


2459398.08 



References. — http://exoplanets.org/almamacframe.html 



-50- 



Table 2. Notations 



Symbol 


meaning 


C 


-5a J (4arf) 


Id 


{aJady/\MjMd) 


e 


1 •^m 1 


V 


Wc - rud 


tc 


{A/3n,){M,/M,){ad/a,f 


X 


ejed 




values of x at centers of librating orbits 



AN drj/dr caused by the disk potential (eq. [A3]) 



